Kitaev model and dimer coverings on the honeycomb lattice 
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We consider an extension of the Kitaev honeycomb model based on arbitrary dimer coverings 
satisfying matching rules. We focus on three different dimer coverings having the smallest unit cells 
for which we calculate the ground-state phase diagram. We also study one- and two- vortex properties 
for these coverings in the Abelian phases and show that vortex- vortex interactions can be attractive 
or repulsive. These qualitative differences are confirmed analytically by high-order perturbative 
expansions around the isolated-dimer limit. Similarities and differences with the original Kitaev 
honeycomb model are discussed. 
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I. INTRODUCTION 

Topologically ordered phases are fascinating states of 
quantum matter characterized by a ground-state degen- 
eracy which depends on the genus of the considered 
topology ' . They are also known to exhibit exotic el- 
ementary excitations, dubbed anyons, which have non- 
trivial braiding statistics"^"^. Recently, the interest for 
the study of topologically ordered phases has been trig- 
gered by the possibility to use anyons as an essential in- 
gredient either to perform quantum computation or to 
provide a quantum memory . To understand these prop- 
erties, much focus has been drawn on two-dimensional 
quantum spin systems especially since Kitaev has pro- 
posed two exactly solvable models, the toric code and 
the honeycomb model , in which such anyons are present. 

In the present work, we provide an extension of the 
Kitaev honeycomb model for which many experimental 
realizations in optical lattices^^^ , Josephson junctions 
arrays^^, or solid state systems with strong spin-orbit 
coupling^ ^'^^ have been proposed. Up to now, many vari- 
ations of this model have already been investigated, but 
almost all of them rely on changing the geometry of the 
underlying lattice^^"^^. Here, we choose an alternative 
route by keeping the honeycomb geometry but by chang- 
ing the coupling configuration. Doing so, we propose a 
set of exactly solvable models which, as we will see, can 
give rise to surprising properties as compared to the orig- 
inal one . 

This paper is organized as follows. In Sec. II we re- 
call the basic properties of the Kitaev honeycomb model 
and we explain how to extend it by considering differ- 
ent dimer coverings. In Sec. Ill, we focus on three simple 
coupling patterns for which we compute the ground-state 
(vortex- free) phase diagrams using Majorana fermioniza- 
tion and analyze the influence of a three-spin interaction 
term. The following Sec. IV is dedicated to sparse vor- 
tex configurations. There, we use a perturbative treat- 
ment to show that the nature (attractive or repulsive) of 
vortex-vortex interaction is very sensitive to the dimer 



configuration. 



II. KITAEV MODEL AND DIMER COVERINGS 

The Kitaev honeycomb model describes a set of spins 
1/2 located at the vertices of a hexagonal lattice which 
interact via the following Hamiltonian 



H 



Yl Yl ^«^f<' 



(1) 



=x,y,z CK — links 



where af are the usual Pauli matrices at site i. In the 
original model , Kitaev considered the case where the 
type of links {a = x, y, or z) depends only on its orien- 
tation as depicted in Fig. 1. 




FIG. 1: Dimer covering / initially considered by Kitaev^ with 
two sites per unit cell. Vectors ni and n2 define the unit 
cell. Red, green, and blue links represent x,y, and z links 
respectively. 

So far, extensions of this model were based on a dif- 
ferent lattice geometry ' ' or a different underlying 
algebra ' . Sticking to a spin- 1/2 model, the exact solv- 
ability is preserved only if the following two constraints 
are satisfied: {i) each site of the system must be trivalent 
and {a) the three links connected to a given site have to 
be of different types (x, y, and z). Thus, the Hamiltonian 
(1) can be defined and solved on any trivalent graph. 



Here, by contrast, we study the Kitaev model (1) in the 
original honeycomb lattice but for different link distribu- 
tions or dimer coverings (we shall use the expressions link 
distributions or dimer coverings as synonyms, so dimer 
coverings do not have the same meaning here as they 
have in the context of quantum dimer models). More 
precisely, we focus on the three simple dimer coverings 
shown in Figs. 1 and 2. These configurations are among 
the simplest periodic decorations one may build on such 
a lattice. Indeed, the honeycomb lattice has two sites 
per unit cell so that any periodic decoration will have 
2n sites per unit cell. Configurations shown in Figs. 1-2 
correspond to n = 1, 2, and 3 respectively, and the corre- 
sponding models will be consequently referred to as /, //, 
and III. In this terminology, the original Kitaev model 
introduced in Ref. is model /. 





FIG. 2: Left: dimer covering // with four sites per unit cell. 
Right: dimer covering /// with six sites per unit cell. Con- 
ventions are the same as in Fig. 1 and labels (a), (b), and (c) 
refer to the different kinds of plaquettes whose corresponding 
operators are given in Eqs. (4-8). 



The three dimer coverings considered here differ by the 
number and the type of plaquette operators Wp. More 
precisely, for covering / there is only one type of plaquette 
operator 



Ti/-^ ^z ^x ^y ^z ^x ^y 



(3) 



where 1 to 6 are the six spins numbered clockwise around 
any plaquette p, spin number 1 being at the bottom of 
the plaquette. Covering // has two types (a and b) of 
plaquette operators (see Fig. 2) 
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Covering /// has three types (a, b and c) of plaquette 
operators (see Fig. 2) 
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One can readily see that the structure of these opera- 
tors is similar for coverings / and // since they contain 
two Pauli matrices of each kind (x, y and z) whereas, for 
covering ///, each plaquette operator is made of identical 
Pauli matrices. 



III. PHASE DIAGRAMS 



As shown below, the phase diagram strongly depends 
on the dimer covering and displays very different prop- 
erties when switching on a three-spin interaction term 
already introduced by Kitaev' to break the time-reversal 
symmetry explicitely. Such a term is built in the follow- 
ing way: if two dimers (i, j) of type a and (j, k) of type 
P share a common site j, then the three-spin term is de- 
fined by crfcF^-cr^ where 7 is the direction orthogonal to 



a and p. 
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Here, we only consider positive couplings J^ and K. A 
crucial feature of Hamiltonian (2) is that it commutes 
with plaquette operators defined, on each plaquette p, as 



Wr, 



a^; 



out(i) 



where i runs over the set of six spins 



around the plaquette p, and where the notation out(z) 
means the "outgoing" direction at site i, with respect 
to the plaquette's contour- . Since Wp = I, plaquette 
operators Wp have eigenvalues Wp = ±1. Consequently, 
there is one conserved Z2 degree of freedom per plaquette. 
Following Kitaev^, we will say that there is a vortex on 
a plaquette if Wp = —1 and no vortex if Wp = +1. Since 
[H, Wp] = 0, one can block-diagonalize H in each vortex 
sector given by a configuration of the Wp^s. 



For the three dimer coverings introduced in Sec. II, we 
shall now determine the ground-state phase diagram. As 
shown by Lieb in a different context , the ground state, 
at least for K = 0, lies in the vortex- free sector {wp=-\-l 
for all p) on which we focus thereafter. 



A. Majorana fermionization 

Undoubtedly, one of the most remarkable properties 
of the Hamiltonian (2) is that, in each vortex sector, it 
can be mapped onto a free Majorana fermion problem in 
a static Z2 gauge field. As detailed in Ref. , H can be 
written in terms of Majorana fermion operators Cj at site 
j, and reads 



H — 7 ^AjkCjCk, 



(9) 



j,k 



where Ajj. = 2Jjj^Ujk + 2KujiUki- In the latter expres- 
sion, Jjk is the coupling on the link (j, k) and Ujk = 
—Ukj = ±1 is a gauge field which describes the vor- 
tex configuration through the relation Wp = Y\(j k)ep^3k 
where j belongs to the black sublattice and k to the white 
one (see Figs. 1 and 2). In Eq. (9), the sum runs over 
all sites and Uij is only nonvanishing when i and j are 
nearest neighbors. 



Thus, for a given vortex configuration, one simply has 
to diagonalize the matrix A to get the spectrum of the 
corresponding sector. Obviously, such a task is easy for 
periodic configurations but difficult otherwise. In partic- 
ular, when one is interested in sparse vortex configura- 
tions, one needs to diagonalize large systems and perform 
some finite-size analysis (see Sec. IV). 

To determine the phase diagram, we computed the 
fermionic gap for the three dimer coverings /, //, and ///. 
For symmetry reasons, all results below are displayed in 



the plane Jx -\- Jy -\- Jz = ^ which is further para 


imetrized 


in the following way 




X = I- Jx^ Jy, 


(10) 


y ^ ^ 'Jx 'Jy ^ 'Jzi 


(11) 



with X, y G {0, 1}. 




B. Dimer covering / 

The dimer covering shown in Fig. 1 is completely 
symmetric under the exchange of the three directions 
x,?/, and z. Thus, for convenience and without loss of 
generality, let us assume J^ ^ Jx^Jy The fermionic gap 
as a function of K evolves as follows. 

• For K = 0, one has 

/^i{K = 0) = Max{2( J, -Jx- Jy), 0}, (12) 

so that one has to distinguish between a gapped phase 
A for Jz > Jx^Jy and a gapless phase B for Jz < Jx^Jy 

• In the small K limit, a gap opens up in the B phase 
whereas, in the A phase, the gap does not depend on 
K. For all K^ the spectrum remains gapless on the line 

J z ^ J X ~r Jy 

• In the large K limit, the gap in the B phase reaches 
the following asymptotic value 



/\l{K = 00)=2\Jz-Jx-JyV 



(13) 



The evolution of the phase diagram when varying K is 
given in Fig. 3. 

As shown by Kitaev, phase A is Abelian, whereas phase 
B contains non- Abelian excitations in the presence of a 
three-spin term . In both cases, anyonic excitations are 
vortices living on plaquettes. 



FIG. 3: Evolution of the fermionic gap in the vortex- free sec- 
tor for covering / as a function of K. Gapped regions dX K — 
correspond to phases A and the gapless domain is identified 
with phase B. For finite K the phase B becomes gapped and 
the gap saturates in the large K limit. 



the following. 

• For i^ = 0, the phase diagram is exactly the same 
as for dimer covering /. Since model // is equivalent 
to model / when J^ = Jy = or Jx = Jy = Jz^ one 
can conclude that the statistics of the excitations in 
phases A and B will also correspond to the statistics for 
covering /. The same argument will hold for model ///. 

• For i^ 7^ 0, the situation is more complex and one 
must distinguish between two regions. 

- For Jz > Jx ^ Jyi the gap is independent of K and 
is still given by Eq. (12). 



For Jz < Jx 



Jy^ the gap has a nontrivial in- 



dependence. It vanishes on the lines 
Jl + 16/^2 ^ ( j^ _ j^ 



(14) 



• In the large K limit, the gap for Jz > Jx^Jy remains 
unchanged whereas the gap in the B phase is made up 
of two planes which intersects on the line Jx ^ Jy = 2/3 
(see Fig. 4). 



C. Dimer covering // 

Contrary to covering /, covering // is not symmetric 
under the exchange of the three directions x, ?/, and z, 
but it is only symmetric under the exchange x -n^ y. In 
this case, the evolution of the fermionic gap with K is 



D. Dimer covering /// 

This covering is completely symmetric under the ex- 
change of the three directions x, y, and z. The fermionic 
gap as a function of K is however very different from 
previous cases, and is illustrated in Fig. 5. 





FIG. 4: Evolution of the fermionic gap in the vortex- free sec- 
tor for covering // as a function of K. For K — 0^ the gap is 
the same as for covering /. 



• For K = 0, the gap is given by 



A///(K = 0) = 2^ J2 + J2 + J2 _ J^Jy _ JyJ^ _ J^J^, 

(15) 
As can be easily checked, it only vanishes at the isotropic 
point Jx = Jy = Jz (as besides for any possible covering 
of the honeycomb lattice at i^ = 0). In other words, for 
this covering there is no gapless phase but just a gapless 
point. 

• For i^ 7^ 0, a gapped non-Abelian phase starts to 
develop around the isotropic point delimited by a gapless 
circle. The gap in the Abelian region decreases when 
increasing K. 

- For K = >/3/18, the gapless circle reaches the middle 
of the triangle sides and isolates each Abelian gapped 
phase. 

- For K = a/3/9, the gapless line finally reaches the 
triangle corners so that only the non-Abelian phase 
remains. 

• In the large- i^ limit, the gap becomes constant and 
equal to Ajjj{K = oo) = 2/3 in the whole parame- 
ter range considered and the excitations are non-Abelian 
anyons. 



IV. VORTEX PROPERTIES 

Let us now discuss the one- and two vortex sectors. 
Our main motivation for such a study is to unveil some 
qualitative differences of the vortex properties for the dif- 
ferent coverings. Such differences are already present in 



FIG. 5: Evolution of the fermionic gap in the vortex- free sec- 
tor for covering/// as a function of K. For /T = 0, there is 
a single gapless point and the whole system is in an Abelian 
phase. By contrast, in the large K limit, the system is in a 
non-Abelian phase and the gap is constant. 



the Abelian phase {J^ > Jx -^ Jy) and for /f = which 
is the parameter range to which we restrict in the fol- 
lowing. Furthermore, for the sake of simplicity, we set 
Jx — Jy = J and Jz = 1 — 2 J. Along this line, vortex 
properties of coverings / and // are identical but differ 
from those of covering ///. 



A. Exact diagonalization 

For this study, we have also used Major ana fermion- 
ization, but contrary to the vortex-free sector, we now 
have to deal with non-periodic configurations and to per- 
form some finite-size analysis on large systems. Here, we 
treated system sizes up to A/" = 16200 x 16200. 

One important remark is that, using periodic bound- 
ary conditions, vortices can only be created in pairs . 
Thus, to compute the one- vortex gap, one has first to 
consider the two- vortex gap A£'2v"' ^^ a pair of vortices 



on plaquettes i and j defined by 



8,. 



(16) 



The right-hand side is the energy difference between the 
ground-state energy £^(*'-^) of the considered two- vortex 
sector and £q which is the ground-state energy of the 
vortex- free sector. Then we use the fact that, in the 
Abelian phase, vortex excitations are gapped so that 
the interaction of two distant vortices decreases expo- 
nentially. This allows us to define the one-vortex gap 
as 



A^iv 



/^E. 



2v 



for \i-j\ >1, 



(17) 





FIG. 6: Effective lattices obtained in the isolated- dimer limit 
J <^ Jz hy replacing z-dimers by effective sites (filled black 
circles). The effective lattice for coverings / (left) and // 
(middle) is a "square" lattice. For covering /// (right), one 
gets a Kagome lattice. Gray-shaded areas show the different 
types of effective plaquette identified by the same labels as in 
Fig. 2. Red (green) bonds denote x (y) links. 



where it is essential that both vortices are of the same 
type, i. e. a or b for model //, and a or b or c for 
model ///, see Sec. II). Practically, results are converged 
for \i— j\ — 10 as long as the vortex gap is not too small, 
namely, away from the transition points (see Sec. IV C). 



B. Pert ur bat ive approach 

To have a better understanding of these few- vortex 
properties, we have also developped an alternative ap- 
proach. We computed the spectrum of H perturba- 
tively up to high orders using the perturbative continuous 
unitary transformation (PCUT) method around the 

isolated-dimer limit, where one of the couplings is much 
larger than the others. Such a strategy has been shown 
to be very successful in the original Kitaev model'^'^''^^'^^ 
although restricted to Abelian phases. Here, we proceed 
along the same line as in Ref. and we skip all technical 
details which can be found in this reference. 

As discussed in Refs. ,23,31,3^, in the limit J <^ J^ 
it is convenient to analyze the problem by replacing z- 
dimers by effective sites with four degrees of freedom 
corresponding to the four possible dimer states. Then, 
depending on the dimer configuration, one can face dif- 
ferent structures (see Fig. 6). 

For covering /, the effective lattice is a "square" lat- 
tice in which all plaquettes are identical. For covering //, 
one still has a square lattice but one must distinguish be- 
tween two kinds of plaquettes (a and b) which reflect the 
operator structure (4). By contrast, the effective lattice 
for covering III is a Kagome lattice with three different 
types (a, b and c) of plaquettes: one hexagonal plaquette 
with alternating x- and ?/-links and two triangular pla- 
quettes with either only x- or only ?/-links. It is worth 
noting that such an effective lattice was also obtained in 
Ref. 17, starting from a different initial geometry. 

As already explained, in the Abelian phase the Kitaev 
honeycomb model contains two different kinds of exci- 
tations: low-energy vortices and high-energy fermions. 
The PCUT method around the isolated-dimer limit pro- 



vides (order by order in J/Jz and in the thermodynamic 
limit) an effective Hamiltonian which commutes with the 
number of fermions. Thus, the low-energy spectrum is 
obtained by considering the zero-fermion (quasiparticle) 
block usually denoted OQP. In this model, the effective 
Hamiltonian in the OQP sector is particularly simple be- 
cause it can be fully expressed in terms of the conserved 
Z2 plaquette operators Wp. Its general structure reads 



H 



OQP 
eff 



En 



Cti 



■ ,P„Wp, 



■ ■W„ 



(18) 



{Pl,--jPn} 



where ^0 is a constant and where {^1,^2, • • • ,Pn} denotes 
a set of n plaquettes. The sum runs over all possible 
plaquette numbers and configurations. In the present 
work, we have computed the corresponding coefficients 
^Pi,...,Pn of ^^is multi-plaquette expansion up to order 
10. The ground-state energy of any vortex configuration 
is then readily obtained by replacing Wp by —1 if there 
is a vortex and by +1 otherwise. 

Before presenting the results, we would like to point 
out that the validity range of the PCUT results does not 
depend on the vortex filling. However, this perturbative 
treatment is only valid as long as all fermionic gaps are 
finite. For instance, this implies that although the gap- 
less phase is reduced to a single point for dimer covering 
1/7, this perturbative approach breaks down well before 
this point is reached since some fermionic gaps in other 
vortex sectors vanish . For the three dimer coverings 
considered here, we are led to conjecture that there al- 
ways exists a vortex configuration such that the fermionic 
gap vanishes for J^ = Jx^Jy For covering / and //, this 
is achieved in the vortex- free sector but, for covering ///, 
this is found for more complex vortex configurations^'^. 
As a consequence, we restrict our analysis to J G [0, 1/4]. 



C. Single-vortex properties 

On the line considered here {Jx = Jy = J), cover- 
ings / and // are identical. However, for covering 1/7, 
one must distinguish between two different vortex gaps 
according to the kind of plaquette which is excited (tri- 
angle or hexagon of the effective Kagome lattice). The 
single-vortex gaps obtained by Majorana fermionization 
and by the PCUT are displayed for these three different 
cases in Fig. 7. An excellent agreement between both 
approaches is evidenced. 

The most striking result is that the value of this vortex 
gap depends strongly on the type of effective plaquette. 
This can be understood easily by noting that the leading 
order in the perturbative expansion is determined by the 
"perimeter" of the plaquette (3 for the triangle, 4 for the 
square, and 6 for the hexagon). In this sense, one may 
say that the mass of the vortex is strongly sensitive to 
the dimerization pattern. 

For illustration, we give the analytical expression of 
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Finally, let us remind that in this Abelian phase, vor- 
tices behave as Abelian anyons or hard-core bosons de- 
pending on their relative position (see Ref. for a de- 
tailed discussion). In the present case, for covering //, 
plaquette excitations (a) and (b) are mutual semions 
whereas they behave, individually, as hard-core bosons. 
For covering ///, triangular and hexagonal plaquette ex- 
citations are mutual semions whereas triangular (hexag- 
onal) plaquette excitations behave individually as hard- 
core bosons. 



D. Vortex- vortex interactions 

To conclude this section, let us discuss the nature of 
vortex- vortex interactions 'y^*'-^) between two vortices on 
plaquettes i and j. The interaction energy between two 
such vortices is defined as 



/i,i) ^ ^^; 



ihj) 



AE^^ 



A^it\ 



(22) 



where the energy of each type of vortex is subtracted. 

Of course, such an interaction depends on the relative 
position of the vortices. Once again, on the line con- 
sidered here {Jx = Jy) this interaction is the same for 
coverings / and //. For these two coverings, we found 
that this interaction is always attractive (this is also true 
for Jx 7^ Jy) whereas, for covering ///, there are some 
positions for which it is repulsive. These results are il- 
lustrated qualitatively in Fig. 8. 



FIG. 7: Single- vortex gap AEi^ as a function of J for cover- 
ings / and // (top) and for covering /// (middle and bottom). 
We give the bare series expansion computed up to order 10, 
some Fade approximants of these series (except for the hexag- 
onal plaquette for which we do not get enough terms in the 
series to provide reliable approximants), and the exact diag- 
onalization results. Insets show the excited plaquette consid- 
ered together with the effective lattice. 



these various gaps computed up to order 10 
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FIG. 8: Schematic view of vortex- vortex interactions v^^'-^^ for 
coverings / and // (left) and for covering /// (right). A green 
(red) plaquette interacts attractively (repulsively) with a vor- 
tex placed on the central blue plaquette. The color code is 
such that a darker green corresponds to a stronger attractive 
interaction. 

In addition. Fig. 9 shows a comparison between ex- 
act diagonalizations and PCUT results for some special 
cases where vortices are placed on nearest-neighbor effec- 
tive plaquettes. Again, both approaches are in very good 
agreement and confirm that some repulsion can be found 
in covering ///. This interesting finding can be readily ex- 
tracted from the series expansion obtained from PCUT. 
Indeed, the leading order of the interaction between two 
effective triangular plaquettes is given by the repulsive 
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FIG. 9: Vortex- vortex interaction v^'^'-^^ for two adjacent effec- 
tive plaquettes as a function of J for coverings / and // (top) 
and for covering /// (middle and bottom). As in Fig. 7, we 
do not show Fade approximants when hexagonal plaquettes 
are involved since the corresponding series starts at order 6 
so that there is not enough terms to get reliable results when 
resumming. Insets show the nearest-neighbor excited plaque- 
ttes for which the two- vortex interaction is computed. 



teraction between a triangular and a hexagonal plaquette 
is attractive and reads — 7J^/256. 



V. SUMMARY AND PERSPECTIVES 

We have extended the Kitaev honeycomb model to ar- 
bitrary dimer coverings satisfying matching rules. This 
results in a family of models which can still be solved ex- 
actly by IVtajorana fermionization. Although these mod- 
els are equivalent at the isotropic point, they strongly 
differ at any other point in parameter space. 

In the present work, we focused on three different 
coverings (/, //, and ///) having the smallest possible 
unit cells. For these three models, we have studied the 
ground-state phase diagram as a function of the three- 
spin interaction strength K. We found that there are 
many similarities between the phase diagrams of Hj and 
Hii. However, the phase diagram of Hm is dramatically 
different. In particular, the gapless phase at K = is 
reduced to a single point and, for large K, no Abelian 
phase exists. We have furthermore studied one- and 
two-vortex properties in these different coverings in the 
Abelian phase at i^ = and we have shown that they 
also depend on the covering considered. These findings 
have been analytically confirmed by perturbative expan- 
sions around the isolated-dimer limit inside the Abelian 
phase. 

Finally, let us point out two routes left for the future. 
Firstly, it would be interesting to investigate other vor- 
tex congurations since, as recently observed in the orig- 
inal Kitaev model , vortex interactions could lead to 
new topological phases characterized, notably, by exotic 
Chern numbers^''^'^''^^'^^. Secondly, it would be worth 
studying the effect of disordered dimer and/or vortex 
configurations to see if an Anderson-type transition can 
occur in such systems. Note that other kinds of disorder 
have recently been considered in Refs. , . We hope 
that the present work will stimulate further investigations 
in a near future. 
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